Effect of coexisting order of various form and wave vector on low-temperature 
thermal conductivity in d-wave superconductors 
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In light of recent experimental evidence of density wave order in the cuprates, we consider a 
phenomenological model of a d-wave superconductor with coexisting charge, spin or pair density 
wave order of various form and wave vector. We study the evolution of the nodal structure of the 
quasiparticle energy spectrum as a function of the amplitude of the coexisting order and perform 
diagrammatic linear response calculations of the low-temperature (universal-limit) thermal conduc- 
tivity. The work described herein expands upon our past studies, which focused on a particular 
unit-cell-doubling charge density wave, generalizing our techniques to a wider class of coexisting 
order. We find that the question of whether the nodes of the d-wave superconductor survive amidst 
a reasonable level of coexisting order is sensitive to the form and wave vector of the order. However, 
in cases where the nodes do become gapped, we identify a signature of the approach to this nodal 
transition, in the low-temperature thermal conductivity, that appears to be quite general. The 
amplitude of this signature is found to be disorder-dependent, which suggests a connection between 
the presence of coexisting order in the underdoped cuprates and recent observations of deviations 
from universal (disorder-independent) thermal conductivity in the underdoped regime. 

PACS numbers: 74.72-h, 74.25. Fy 
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I. INTRODUCTION 

The low energy excitations of cuprate superconduc- 
tors are Dirac fcrmions, which arise due to the d-wave 
nature of the superconducting order parameter. One 
expected signature of these nodal quasiparticles is the 
presence of a universal term in the low-temperature ther- 
mal conductivity, koo, which depends only on the ratio 
of the gradient of quasiparticle dispersion to the gra- 
dient of the gap, Vf/vA, but not on the disorder^! 3 .. 
In the optimally doped and overdoped regimes, koo has 
been observed in many instances, and agrees closely with 
its predicted V aluei«M»L2»2»12»ii 1 i2 i 12»l£. For instance, in 
YBa2(Cui_ x Zn x )306.9, koo is observed to be insensitive 
to the concentration of Zn atoms, which are varied to 
allow the scattering rate to range over several orders of 
magnitude 4 . However, in some cases, the value of the uni- 
versal limit thermal conductivity, kqq, is different than 
expected, or not observed at all, in particular as one 
approaches the underdoped regim o 15 ' 16 i 17 i 18 i 19 i 21 . One 
possible reason may be that disorder-induced local mag- 
netic moments enhance the scattering rate while leav- 
ing the density of states unaffected, thus reducing the 
transport^. 

Another mechanism which might account for devia- 
tions from the expected value of kqq is the presence 
of competing order parameters. For years, evidence 
of the presence of additional symmetry breaking or- 
der parameters in cuprates has been compiled in neu- 
tron scattering data and scanning tunneling microscopy 
c::pcrimcnt a 22 - 23 - 24 - 25 - 26 - 27 - 28 - 29,30 - 31,32 - 33,34 - 35 - 36 - 37 The 
presence of additional orders other than superconduct- 
ing may be incidental, yet it also may be intrinsically 
related to the complex phenomenon of high temperature 
superconductivity itself 38 . 



The addition of order parameters which preserve time- 
reversal symmetry followed by a lattice translation was 
found to preserve the nodal structure of the quasipar- 
ticles, for small amplitudes of order—. As the strength 
of such ordering perturbations increases, the locations of 
the nodal excitations evolve in fc-space. For sufficiently 
large amplitude, the quasiparticle spectrum can be en- 
tirely gappe d 39 ' 40 ' 41 ' 42 . Such a modification of the quasi- 
particle spectrum should manifest itself in the low tem- 
perature thermal conductivit y 41 ' 42 ' 43 . 

In this paper, we model a cuprate superconductor us- 
ing a mean-field formalism describing a BCS-like e?-wave 
superconductor (dSC) perturbed by an additional order 
parameter. We calculate the low-temperature thermal 
conductivity, accounting for the presence of several dif- 
ferent varieties of competing orders. We argue that these 
predictions can then be used as an indirect verification 
of the presence or absence of various competing orders in 
cuprates. 

A previous linear response calculation of kqq in a dSC 
with the addition of a Q = (tt, 0) charge density wave 
showed that vertex corrections were not important for 
the universal limit thermal conductivitity, within the self- 
consistent Born approximation 4 ^. As the charge density 
wave's amplitude increased beyond a critical value ip c , 
the quasiparticle spectrum became gapped. Correspond- 
ingly, the thermal conductivity (made anisotropic by the 
presence of the density wave) vanished beyond that crit- 
ical strength of ordering. In addition, a dependence on 
disorder resulted, in particular for charge orderings near 
the transition. Armed with this information, we proceed 
to study the effects of a wider class of density waves on 
the low-energy properties of cuprates. 

In Sec. [Ill we will develop the mean-field formalism 
we will use to describe superconductors in the presence 
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of competing orders. We write effective hamiltonians 
for charge, spin and pair density waves of several wave- 
vectors. Additionally, we describe configurations with 
multiple wave vectors, such as checkerboard order. Next, 
in Sec. lIIII we will derive the current operators associated 
with the various kinds of orders, and use this to estab- 
lish a relation for the bare-bubble thermal conductivity. 
Finally in Sec. I1V1 we will apply our results to several 
different cases, and compare and contrast the results. 



Upon Fourier transform this becomes 

H C dw = ^2(®Qfkc{ +Qa c ka + $* Q fkc{a c k+Qv), (7) 



with the definition $q = ipe 1 ®'*", where r*o describes the 
shift of the density wave from being site-centered and tp is 
the amplitude of the density wave. Coupled with Eq. ([5]), 
this indicates restrictions on certain density waves' reg- 
istration with the lattice. 



II. MODEL 

A. States of broken symmetry 

States that arise as a result of broken symmetries 
are characterized by the presence of non-vanishing off- 
diagonal matrix elements. The superconducting state it- 
self, for instance, can be identified with the non- vanishing 
anomalous Green functions, as was shown by Gor'kov 44 . 
These anomalous Green functions are defined in space 
and time as 



< VaM)#(0,0) > 



(1) 



Given singlet paired electrons of opposite momenta, this 
corresponds, in momentum space, to 



< Va(fc,t)#(-fc,i) > 



(2) 



In a similar fashion, ordered states representing density 
waves will also admit non- vanishing correlations between 
states separated by the wave- vector of the density wave. 
In this chapter, we will consider the subset of those which 
are defined in momentum space as 

< V Qt (fc + Q,i)V/3(M) >= $> Q f(k)d(a,f3), (3) 

representing charge (d — 5^) and spin (d = 5p{5" — 8")) 
density waves, as well as pair density waves 



< 



(4) 



For the purposes of simpler classification of orders, we 
are carefully following some definitions made by Nayak 
in Ref. |45|, so that &q will represent the magnitude and 
phase of the density wave and f(k) is an element of a 
representation of the space group of Q on a square lattice. 
Certain order parameters obey restrictions. For instance, 
charge and spin density waves for which 2Q is a member 
of the reciprocal lattice obey the additional condition 



(5) 



as was pointed out in Ref. |45| . 

Written as a sum over real space, the hamiltonian rep- 
resenting a charge density wave system is 



Hcdw = V>e- i<5 ' (F - Fo) /(r - ri)c\ a c Tla + h. 



(6) 



B. d-wave superconductor 

Our starting point is a model for d-wave superconduc- 
tors 



H = J2 { e kcl a c ka + A fc 4 T cL fcl ) + h.c. 



(8) 



where the normal state dispersion is given by a tight- 
binding hamiltonian, 

£k = — 2i(cos k x + cos k y ) — t' cos k x cos k y — [i. (9) 

and the superconducting order parameter is of d x 2 -y 2 
symmetry, 



A fc = ^(cos(M - cos(M). 



(10) 



As given, this hamiltonian has nodal excitations, which 
are located along the d x 2_ y 2 symmetry axis in the 
(±7r, ±7r) directions. The nodes' distance from the 
(±7r/2, ±7r/2) points is controlled by the chemical poten- 
tial \i. These quasiparticles are massless Dirac fermions 
in the sense that they have conical dispersion. The exci- 
tation energy is 



E k 



A 2 



(11) 



and at low energies, ~ f/fci and Aj. ~ v&k 2 , where 
k\ and k 2 are momentum-space displacements from the 
node in directions perpendicular and parallel to the Fermi 
surface respectively, Vf is the Fermi velocity, and «a is 
the slope of Afc at the node. For n on the order of t 
or smaller, the ratio of Fermi velocity to gap velocity is 
given as 



t 2 - ft' 2 



A 



(12) 



Then, as perturbations are turned on, the locations of 
the nodes evolve in fc-space, while the stability of the 
nodes is generally preserved for non-nesting perturba- 
tions which preserve the composite symmetry of lattice 
translation followed by time-reversal, as was noted by 
Berg and Kivelson^. 
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C. Density waves of different wave vectors 

The presence of a uniform density wave in a supercon- 
ductor changes the system in both real and momentum 
space. In real space, the unit cell increases. In momen- 
tum space, we see an effective reduction of the Brillouin 
zone, also called band folding. Accordingly, our sec- 
ond quantized descriptions of the systems are modified. 
Whereas in a superconductor we can rewrite a quadratic 
effective hamiltonian using Nambu field operators, 



4 = (4 



i'k = ( j 
■ c -ki 



(13) 



we can alternatively write extended Nambu vectors, such 
as 



4 = (4 



fcf C-kl C fc+Q| C-k-Ql 



(14) 



where the wave- vector Q connects each point in the first 
reduced Brillouin zone with a point in the second re- 
duced Brillouin zone. Sums over fc-space are then per- 
formed by integrating over the reduced Brillouin zone 
(the shaded regions in FigfTJ) , and taking the trace of the 
now-extended Nambu space matrix. The two descrip- 
tions are equivalent, but the extended Nambu descrip- 
tion naturally fits the effective hamiltonians of systems 
with non-zero mean-field density waves. In Fig. [H we il- 
lustrate four different density waves which are considered 
in this paper: Q = (tt, 0), Q = (7r/2,0), Q = (n 7 Tr) and 
Qi = (it, 0), Q2 = (0, 7r) (checkerboard) orders. These 
disturbances are illustrated in real space in Figf2] 



1. Q = (tt, 0) density waves 

A density wave of wave vector Q = (tt, 0) corresponds 
to a striped system: the unit cell is doubled in the x- 
direction, and the Brillouin zone is reduced by 50% as 
seen in Fig. [Ha). The extended Nambu vector is that of 
Eq.[H with Q = (?r,0). 



2. Q = (tt/2,0) density waves 

A density wave of wave vector Q = (tt/2,0) corre- 
sponds again to a striped system, one in which the unit 
cell has increased by a factor of four, and the Brillouin 
zone is reduced by the same factor. The reduced Brillouin 
zone is taken to be the region containing the pre-density- 
wave nodal quasiparticle excitations of the d-wave super- 
conductor; in Fig. [T](b) this is indicated with shading. 
The extended Nambu vector is that of Eq. IT71 



3. Qi = (tt,0): Q2 = (0, tt) checkerboard density waves 

Two density waves of equal weight in orthogonal di- 
rections corresponds to a checkerboard ordered system. 



a) 



b) 



d) 



FIG. 1: Illustrated are the reduced Brillouin zones of a square 
lattice in the presence of different density waves. The dots il- 
lustrate approximately the location of nodal excitations in the 
original d :r 2_ !/ 2-symmetry superconductor; the dashed line is 
the new zone boundary induced by the density wave. Illus- 
trated are density waves of wave vector: (a) Q = (tt,0), (b) 
Q = (tt/2,0), (c) Qx = (tt,0), Q 2 = (0,tt) and (d) Q = (tt, tt) 



As in the case of the Q = (tt/2,0) case, the Brillouin 
zone's area is reduced by a factor of four, although it is 
a different reduced Brillouin zone, illustrated in fc-space 
in Fig.QJc), and in real space in Fig.^c). 



4- Q = (tt,Tt) density waves 

A density wave of wave vector Q = (tt, tt) corresponds 
to a system which is modulated in both the k x and k y 
directions: Sp oc sin(k x ) sin(k y ). The reduced Brillouin 
zone is indicated in Fig. QJd) as the shaded region, and 
the real space modulation is illustrated in Figf5)(d). The 
extended Nambu vector is as in Eq. 1 141 with Q now rep- 
resenting the (tt, 7r) density wave. 



D. Charge density waves 

A commensurate charge density wave is one for which 
the charge density is oscillatory in real space and repeats 
itself after translation by an integer number of lattice 
constants. The momentum space description of a mean 
field hamiltonian for such a system is 



CDW 



ha 



t 

k+C 



T Cka 



(15) 
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b) 5p~Cos(7t/2 x) 
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A charge density wave p = po + Sp which doubles the unit 
cell (so that Sp alternates sign from cell to cell in the x- 
direction) has wave- vector Q = (71", 0). A (tt, 0) CDW 
perturbation in its 4-componcnt extended Nambu basis 
(particle, hole, shifted particle, shifted hole) is given by 



H C dw = 









4* 














-A 


A k 























V-fc, (16) 



c) 8p~Cos(ii x)+Cos(7t y) d) 8p~Cos(7i x)Cos(7i y) 



FIG. 2: Illustrated are the four different density waves con- 
sidered in this chapter, in real space. Each circle corresponds 
to the position of a Cu atom, and the size of the circle indi- 
cates whether the density at that site is higher or lower than 
the average. Illustrated are density waves of wave vector: (a) 
Q = (tt.O), (b) Q = (tt/2,0), (c) Qx = faO), Q 2 = (0,tt) 
and (d) Q = (tv, tv) 



where the sum is over the reduced Brillouin zone of 
Fig.EEta), and we define A k = $ Q fk+$* Q fZ + Q- A (^A 0) 
CDW perturbation written in its 8-component extended 
Nambu basis 



is written as 



W; — ( C [-T C ~ k i C l+2QT C -fe-2Qi C fe+Qt C ~ k ~Ql C /c+3QT 



where H k is given by 
/ 



; I+3QT c -k-3Qlj 



H. 



CDW 



A k 





A* 










-A-k-Q 





-AU- 


-2Q 


4* 

^fc+3Q 





A k +2Q 







V 


4* 





-A^ 


-3Q 



A* k A k+3Q 

-A*_ k - Q 

Ak+Q A* k+2Q 
-A_ fe _ 2Q 



\ 


A* 

^-k-3Q 



J 



(17) 



(18) 



(19) 



E. Pair density waves 

Scanning tunneling microscopy experiments have re- 
vealed the presence of modulations in the local den- 



sity of states in the vortex cores of the cuprate 
Bi2Sr 2 CaiCu2 8+ ^^^i21^^, and in some instances, 
in the absence of magnetic fiel d 24 ' 46 . More recent mea- 
surements, conducted in the absence of magnetic field, 
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measured the spatial dependence of the superconducting 
Their finding was that the superconducting order 
parameter is modulated, corresponding to superconduct- 
ing pairs with a net center-of-mass momentum, which 
has become known as a pair density wave. 

In addition to considering the effects of spatial modula- 
tions in the charge density, we therefore can also consider 
adding (to the dSC hamiltonian) a term corresponding to 
a modulation in the pair density. A pair density wave of 
wave vector Q is written as 

Hpdw = ( Q Q9kcl+Q a c -kp + %9k c -kpc k+Qa y 

k 

a/3 

(20) 

or, with the definition B k = Oq (gk + 9-k~o) we can 
write (for (-7T, 0) or (-7T, tt) density waves) 



Hpdw = 



£4 











B- k 






B*_, 








0/ 



i>k (2i) 



F. Spin density waves 

The effective hamiltonian corresponding to a spin den- 
sity wave of wave vector Q is 



H S dw = ^vy^QfkC 



T Ck0 



•*0/*cLrPk+Q<r 1.(22) 



In the extended Nambu basis, for a Q = (tt, 0) or Q 
(tt, tt) SDW, this takes the form 



H. 



SDW 









A% 











A k 











A* 

A -k 






fa, (23) 



where A k = $ Q fk + %fk +Q - 



G. Checkerboard density waves 



In addition to broken symmetry states arising due to a 
single density wave, we can also consider multiple density 
waves. Scanning tunneling microscopy experiments have 
previously revealed the presence of checkerboard order in 
BiSCCoQ 22 i 23 i 24 i 25 i 26 i 27 . While the wave vectors of the 
order in those experiments was seen to be near Q w 5, 
for simplicity we first write down the hamiltonian corre- 
sponding to Qi = (tt,0) and Q2 = (0, tt) checkerboard 
order. The Brillouin zone is reduced to one fourth of its 
size, as is seen in Fig. [He). The extended Nambu vector 
which describes such a system is 



V4 = (4t c -kl c -k-QA c 1+q bT c -k-Q y l cl+ Ql +Q v ] c_ fc -Q x -Q y |) (24) 

and the second quantized hamiltonian which describes the addition of a charge density wave and pair density wave is 



rrchcckcrboard 
fl CDW 



rrchcckcrboard 
fl PDW 



£ ^k H klpk 



(25) 



where H k is given by 



/ 











B^ 








B 1* 


-a^ 


B -k 


-aH 


A[ x) 
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D k 


-aH 




















A (y) 


B (y) 

h Qx 


A (x) 

^k+Qy 


B (X k Q 


V 





n{y)* 


k Q x 


k-\-Qy 


—A^l* Q 






Ay)* 



B 



(y)* 

-k-Q 
,(x)* 
l k+Q v 



B 



(X 1 c 




B 

-A 
B 
-A 






(y) 

k+Q x 



(y) 

-fc-C 
(x) 

k+Qy 
{x) 



(26) 



r 



and where 



.4 

B 
B 



(x) = 

k — 

(y) = 

k — 
(x) = 
k — 

(y) - 



®Qjk x +$*Qjk* + Q x 

^Qyfky + ^qJk + Qy 

©QxC^ + .9-fc x -Q x ) 

®Qy(9ky + Q-ky-Qy) 



represent the amplitudes of the charge density and pair 
density waves in the x and y directions. 



(27) 
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III. THERMAL CONDUCTIVITY 

At low temperatures, the temperature dependent 
phonon contribution to thermal conductivity vanishes as 

a power, Kphonon ~ ^4,5,6,8,9,11,13,15,16,17,18,19 There _ 

fore, the T-linear quasiparticle current can be extracted 
from experimental data by plotting the measured ther- 
mal conductivity divided by temperature as a function 
of T Q_1 . In previous wor k 41 ' 42 , we considered a site- 
centered Q = (tt, 0) charge density wave and calculated 
the thermal conductivity using Green's functions ob- 
tained from the self-consistent Born approximation, and 
incorporated vertex corrections within the ladder approx- 
imation. Because the results of this work indicated that 
vertex corrections can usually be neglected, in what fol- 
lows we will derive the thermal conductivity using the 
bare-bubble correlation function. This will greatly sim- 
plify the calculation, allowing its application to a vari- 
ety of systems, whereby the thermal conductivity can be 
computed numerically given an an effective hamiltonian 
H = Hdsc + Hnwiip) an d an effective scattering rate Tq. 



Now we have a generalized velocity operator in the 
Nambu space, Vk = For instance, for the Q = (tt, 0) 
pair density wave of Eq. [2TJ the velocity operator would 
read 



v(k) 



( Vf,k 


\ Ml 

\ dk 



VA,k 

dA- k 
dk 






dA*_ 



dAk 
dk 





\ 



dk 

Vf,k+Q V A ,k+Q 
VA,k+Q -Vf,k+Qj 



(32) 



For density waves without internal momentum depen- 
dance, or for those where the variation is slight near the 
nodal locations, the velocity operator reduces to the form 
found in Refs. \MM 



v f,k 



v J,k 
VA,k 



VA,k 



(33) 



B. Universal limit thermal conductivity 



A. Current operators 

In order to calculate the thermal conductivity, we first 
need to derive the heat-current associated with the quasi- 
particles. Because heat and spin currents are both pro- 
portional to the quasiparticle current, we can get the 
heat-current by calculating the spin-current, and then 
using the energy measured from the Fermi level as the 
associated charge (instead of the spin). To calculate the 
spin current for any particular hamiltonian, we write the 
density operator in second quantized form, and then use 
Heisenberg equations of motion to find the momentum 
space representation of the current, that is 



lim(q-j s 

g— 

The density operator is 



= [p S g ,H]- 



(28) 



Pq = J2 (4t c *+«T + c_ fci c T _ fe _ gi ) . (29) 
k 

Taking the commutator with the hamiltonians of Eqs. 
(|8"]) , (fT5|) . and (f2"2"|) , using anti-commutation relations, 
and discarding boundary terms, we find 

[p S q, H ] = Yy^k'^'+qa^kHkM 
kk' 



(30) 



for the spin current. The heat current in the Matsubara 
representation is given by 



k dk 



Y 



(31) 



The universal limit thermal conductivity is calculated 
using linear response formalism. The thermal conduc- 
tivity is given in terms of the retarded current-current 
correlation function. 

K(n,T) Im(n Rct (»)) 

— r— = £3. or* — (34) 

We evaluate the correlation function using the Matsubara 
method 4 ^. The bare-bubble correlator, given in terms of 
a spectral representation, is 

IT(iQ) = / dwidw 2 Tr^ ^A(wx)vA(uj 2 )v S{iQ) (35) 



where 



(iui — uji)(iuj + ifl — UI2) 



(36) 



and A(k,uj) is the spectral function. 

It is important to use the correct form of the spectral 
function in Eq. (|35|) to avoid erroneous results, as is noted 
in Ref. EJ. For example, a bond-centered CDW of wave 
vector Q — (tt,0), which looks like 



H, 



CDW 



/0 iip 
-iip 
#000 

\0 -ii/j 



(37) 



leads to a spectral function which is not real, and the 
spectral function is not given by the formula 



A(k,w) 



1 



(38) 
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but rather by 



-1 



A(k,w) = —[G R (k,u;) 



G A {k,u) 



(39) 



The details of the thermal conductivity calculation are 
similar to those of Refs. Q1 12. In general, the self- 
consistent t-matrix approximation can be used to com- 
pute the Green's functions, however, here we use a sim- 
pler, diagonal self-energy as a first approximation. In 
terms of a model hamiltonian H k , and incorporating im- 
purity scattering by assuming an imaginary part of the 
self-energy, £_r(w — ► 0) = — iTq, the universal limit ther- 
mal conductivity is 



k k 2 B ir 2 r„ dH k ~,„, dH k 

urn — = —2- — 



t^o T 
where 



^ Re [Tr [2(0)^1(0) 



dk 



(40) 



G B (k,w) = (w - ff fc + *T(w)J 
G A (k,u) = (lu - H k - ir(w)) \ 



(41) 



IV. EFFECTS ON ENERGY SPECTRUM AND 
THERMAL CONDUCTIVITY 

Here we modify the dSC hamiltonian of Eq. with 
the addition of density waves such as those of Eqs. (fl~5|) . 
(|20| and (|22| , which will be tuned by the real parameter 
ip, the strength of the density wave. This is done to study 
the behavior of the quasiparticle spectrum, and through 
Eq. (|40|). the universal limit thermal conductivity. In 
each of the figures from Fig. [3] to Fig. ??, we present (a) 
The trajectory of the nodes in the region < k x , k y < 
as the density wave is turned on (the starting place (node 
for dSC system) is indicated with a star), (b) (minimum) 
quasiparticle energy as a function of the order parameter 
strength ip, and (c) universal limit thermal conductivity 
as a function of ip. In all instances, the universal limit 



conductivity is given in units of 



the value 



Zh VfVA 

for the original dSC system, and we measure Ao, E, 
fi and Tq in units of t, the hopping parameter. 



A. Q = (ir, 0) density waves 

The addition of a Q = (ir, 0) charge density 
wave to a d-wave superconductor has been considered 
befor o 40 ' 41 i 42 . As the perturbation is turned on, the 
nodes' locations evolve along curved paths, until they 
meet the images of the nodes from the second reduced 
Brillouin zone at the collision point (ir/2,ir/2), as seen 
in Figs. [3] and 2] . The effect is the same, regardless of 
whether the density wave is of s-wave ($q = ip, fu — F 



site-centered), p^-wave (^q = iip,fk = sin(fc x a), bond- 
centered) or p y -wa,ve ($q = ip, f k — sin(k y a), site- 
centered) symmetry. The critical value of ip which gaps 
the system is ip c = Vfko, where 



k = V2 



cos 



\2t 



t l ( ± )2 _f_ 
y 2t'' At 



(42) 



is the distance separating the ip = nodal point from 
(7r/2,7r/2) in fc-space. For an s-wave perturbation of 
strength ip representing a charge, pair, or spin density 
wave, the quasiparticle spectrum is 



A - y/A 2 + B~C, 



(43) 



where 



A = 
B = 

C EE 



A 2 fe 



A 2 



2ip 2 (ae k e k +Q + bA k A k+Q ) 



(4 + A^)(6 



2 

k+Q 



), 



(44) 



and (a, b) = (—1,1) for a charge density wave, (1,-1) 
for a pair density wave, and (—1, —1) for a spin density 
wave. 

The resulting thermal conductivity is anisotropic, re- 
flecting the striped nature of the system. The nodes are 
deformed as they approach the collision point, and the 
thermal conductivity K yy perpendicular to the direction 
of the density wave increases at first, before both k xx 
and Kyy vanish for larger amplitudes of density wave, ip. 
The effect of a Q = (it, 0) pair density wave is similar to 
that of the charge density wave: the nodes evolve along 
a curved path until they meet their images in the sec- 
ond reduced Brillouin zone, and the resulting universal 
limit thermal conductivity is the same. The effects of a 
site-centered (tt, 0) pair density wave is shown in Fig. [5] 

A more unusual case is that of the Q = (jr, 0) spin 
density wave. With this perturbation, the nodal points 
evolve directly towards the (7r/2, k y ) line, as seen in FigJHl 
The quasiparticle spectrum then evolves so that there are 
two minima. In other words, the node splits in two, and 
nodes move up and down the (ir/2,k y ) line. The nodes 
are nested by Q, but the spectrum remains gapless, and 
the universal limit thermal conductivity is unaffected. If 
the perturbation is allowed to become extremely large 
{ip » Ao), then the nodes (there are now twice as many, 
each on a reduced Brillouin zone edge) collide with their 
images from the second reduced Brillouin zone, and the 
thermal conductivity then vanishes. The split-off nodes 
collide at different strengths of ip, however, and the spec- 
tral weight disappears in two steps, accordingly, as does 
the thermal conductivity. The fact that the nodal struc- 
ture is preserved in this case, even when the nodes be- 
come nested, runs contrary to the intuition (suggested by 
the converse of the theorem of ReflU) that such nested 
nodes would become gapped. 
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FIG. 3: Effects on spectrum and low temperature transport 
of a Q = (tt, 0) charge density wave. The results are the same 
for site centered (s-wave or p^-wave) and bond centered (p x - 
wave) density waves, in that the nodes evolve along the same 
curved paths toward the (±7r/2, ±vr/2) points where they col- 
lide with their images from the next reduced Brillouin zone. 
As this happens, the nodes are nested and the spectrum is 
gapped. The universal limit thermal conductivity vanishes 
beyond this point. Disorder To broadens the transition. Here 
we take fj. = -0.6, A = 4 and F = 0.02. 
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FIG. 4: Effects on spectrum and low temperature transport 
of a Q = (tt, 0) charge density wave. Depicted are the results 
for (i = — 1, Ao = 0.4 and To = 0.02. These parameters 
describe anisotropic Dirac quasiparticles, with Vj/va 
The anisotropy tends to suppress kqo slightly. 
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B. Q = (tt, 7t) density waves 

Adding a Q = (tt, it) spin density wave was also dis- 
cussed as an example in Ref. [39l In real space, such 
a density wave is modulated as cos(k x ) cos(k y ) , so that 
nodes remain along the (it, tt) direction as the density 
wave is turned on, as is seen in Fig. [7] When the nodes 
reach (tt/2,tt/2), the system is gapped 3 -, and the thermal 
conductivity vanishes. On the other hand, a Q = (tt, tt) 
charge density wave behaves in a similar manner to the 
(tt, 0) spin density wave, in that the nodes do not vanish 
for small perturbations. 

The addition of a Q = (tt, tt) pair density wave 
drives the location of the nodes towards the T point at 
(k x , k y ) = (0, 0), an effect which is also observed via the 
addition of the checkerboard pair density wave. In both 



FIG. 5: Effects on spectrum and low temperature transport 
of a Q = (tt, 0) pair density wave. Depicted are the results for 
/i = —0.6, Ao =4 (vf = va) and Fo = 0.02. As was the case 
for the CDW, the nodes evolve along a curved path towards 
the (n/2,k y ) line. Upon reaching k x — tt/2, the nodes are 
nested, and the spectrum is gapped. For tp larger than the 
critical value tp c , the thermal conductivity vanishes, up to 
disorder broadening. 
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FIG. 6: Effects on spectrum and low temperature transport 
of a Q = (tt, 0) spin density wave. Depicted are the results 
for parameters Ao = 4, fi = —0.6 and To = 0.02. As the 
density wave is turned on, the nodes move in a straight line 
to the (±7r/2, k y ) lines. When they reach that line, each node 
splits in two, and the two nodes move up and down along 
that line. The spectrum remains gapless, even though the 
nodes are nested by the ordering vector. Correspondingly, 
the thermal conductivity is unaffected at that energy scale. 
For ip much larger, these two nodes collide with their images 
in the second reduced Brillouin zone (at different values of 
tjj), and the thermal conductivity is reduced by one half of 
the pure dSC value after each such collision. The locations of 
the two separate nodal collisions are illustrated by a square 
and a triangle in (a), and the order strength at which they 
appear is given in (b) and (c). 



instances, a large perturbation -0 >> Ao is required to 
affect the thermal conductivity. 



C. Q = (tt/2, 0) charge density wave 

A Q = (tt/2, 0) density wave behaves slightly differ- 
ently from the (tt, 0) case. In this case, the nodes are 
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D. Qi = (7r,0), Q2 = (0, tt) checkerboard density 



Configurations with more than one density wave can 
also be considered in this formalism. In this paper, we 
turn our attention to the checkerboard configuration il- 
lustrated in part (c) of Fig. [1] As we turn on two charge 
density waves of Qi = (71", 0) and Q2 = (0, 7r), with equal 
amplitudes, the nodes are perturbed along the symme- 
try line toward the (tt/2, 7r/2) point, as shown in Fig. [9] 
When the nodes reach the (tt/2, tt/2) point, the spectrum 



FIG. 7: Effects on spectrum and low temperature transport of 
a Q = (7r, tt) spin density wave. As the density wave is turned 
on, the nodes move along the symmetry lines k x = ±k y to- 
ward (±7r/2, ±7r/2), where they become gapped. Accordingly, 
/too vanishes. In (c) the effects of increasing disorder are pre- 
sented. The disorder tends to smear the thermal conductivity 
around the nodal transition; as such, Koo is no longer univer- 
sal. 



driven towards the (rr/A,k v ) line, rather than (rr/2,ky). 
While they would become gapped if they arrived there, 
for realistic parameters t, fi, and Ao, such a density wave 
would dominate the system, that is, ip >> Aq. The evo- 
lution is as seen in Fig. [5] and preserves the nodes for 
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FIG. 9: Effects on spectrum and low temperature transport 
of a Qi = (7r, 0), Q2 = (0,7r) charge density perturbation to 
the dSC system. The nodes move in straight lines toward the 
(±7r/2, ±7t/2) points. The spectrum becomes gapped at that 
point, and kqo vanishes for ip larger than about OAt. 
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becomes gapped, and the thermal conductivity vanishes, 
with a value of ip c about two-thirds of that for the striped 
(tt, 0) CDW. In contrast, the checkerboard pair density 
wave seen in Fig. |TD] evolves the nodes along the same 
symmetry line, but towards the T point (0,0). At that 
point, the spectrum would become gapped, and the uni- 
versal limit thermal conductivity would vanish. However, 
systems which more closely resemble a d-wave supercon- 
ductor than the checkerboard (Ao > ip) will remain gap- 
less, as the nodal evolution would not be driven that 
far — about thirty times the critical value for the striped 
(tt, 0) PDW . 



FIG. 8: Effects on spectrum and low temperature transport 
of a Q = (7r/2, 0) charge density wave. As the density wave is 
turned on, the nodes move in a curved path to the (±7r/4, k y ) 
lines. The spectrum becomes gapped at that point, when 
the node is nested by the ordering vector and the thermal 
conductivity vanishes for ip larger than about 4t. Such a 
system is out of the range of validity of our model, as it would 
be dominated by the charge order, rather than the d-wave 
superconductor. For ip on the order of ip c , Koo retains its dSC 
value. Here, /1 = —0.6, Ao = 4 and To = 0.05 



ip < Ao- As such, the universal limit thermal conductiv- 
ity is not significantly affected by this perturbation. 



V. CONCLUSIONS 

In conclusion, we have written mean-field hamiltoni- 
ans describing a d-wave superconductor perturbed by a 
variety of density waves. We noted the effects of such per- 
turbations on the low energy quasiparticle spectrum, and 
by calculating the universal limit (T — ► 0, f2 — * 0) ther- 
mal conductivity, see the effect that density waves can 
have on the low temperature thermal transport. Whether 
or not the universal limit thermal conductivity is robust 
in the presence of an incipient density wave depends on 
which type of density wave, and which wave vector, is 
added. For instance, in the case of Q = (tt, 0) pair den- 



10 



3ff 
1 



a) cli 


*ckerb 


oard PDW 















A 



A 




FIG. 10: Effects on spectrum and low temperature transport 
of a Qi = (7r, 0), Q2 = (0, it) pair density perturbation to 
the dSC system. The nodes move in straight lines toward the 
origin. The spectrum would become gapped at that point, 
however, the value of ij) required is much larger than the en- 
ergy scale of the superconducting order parameter. Therefore, 
for reasonable strengths of the ordering vector, the thermal 
conductivity is unaffected by this density wave. 



sity waves, the quasiparticle nodes evolve so that their 
/c-space locations move toward (±7r/2, ±7r/2), as they do 
for a CDW of the same wave vector. When they reach 
this point, which is the point at which the density wave 
vector nests the nodes, the spectrum becomes gapped. 
However, for the Q = [it, 0) spin density wave, the nodal 
structure is preserved beyond this ordering strength, de- 
spite the nesting of the nodes. For the Q = (it, it) den- 
sity waves, the effects of SDW and CDW are reversed 
from that of the (71", 0) case; the (it, it) CDW preserves 
nodality beyond the nesting wave vector, while the (it, it) 
SDW is gapped beyond a critical strength. In the case of 
Q = (7r/2,0), the different wave vector drives the nodes 
toward (7r/4, k y ) instead. Given typical tight-binding pa- 
rameters, such a charge density wave will not gap the 
quasiparticle spectrum, and will thus not affect the ther- 
mal conductivity, which remains universal. In the case 
of Qi = (it, 0), Q2 = (0, it) checkerboard charge order, 
whether or not the universal limit thermal conductiv- 



ity is robust depends on which type of density wave is 
present. The CDW checkerboard nodes move toward the 
(±7r/2, ±7t/2) point, and become gapped. However, the 
PDW checkerboard nodes move away from that direction 
and the nodal structure is preserved. 

Because the onset of charge ordering is believed to be 
correlated with underdoping, observations which show 
that the low temperature thermal conductivity differs 
from the universal value predicted in Ref. [U may be due 
to the influence of coexisting orders. There are some 
general features that appear in all of the models consid- 
ered in this paper. (1) In general, the nodal evolution 
is determined more by the wave-vector Q than by the 
chemical potential /1, although /1 will determine the am- 
plitude of density wave which will gap the system. (2) 
The physics still remains nodal in the following sense. 
Whether the density waves considered were of s-wave or 
p-wave symmetry did not have an effect; all that matters 
is the amplitude of the density wave at the node. (3) 
It is interesting to note that the universal limit thermal 
conductivity generally develops a disorder dependence, 
especially near the nodal transition point. The presence 
of density waves are therefore one possible explanation 
of the breakdown of universal limit thermal transport in 
cuprates. (4) In general, there is an increase in the ther- 
mal conductivity near the nodal transitions (for cases 
where there are), which is caused by the deformation of 
the nodes (and the resulting effective change in U//«a) 
as they meet their images in the second reduced Brillouin 
zone. This feature is consistent with the thermal conduc- 
tivity measurements of Proust et al£^ who find a large 
enhancement in thermal conductivity of Bi-2201. 
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